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Abstract. A finite element method for solving the resonance line transfer problem in moving media is presented. The algorithm 
works in three spatial dimensions on unstructured grids which are adaptively refined by means of an a posteriori error indicator. 
Frequency discretization is implemented via a first-order Euler scheme. We discuss the resulting matrix structure for coherent 
isotropic scattering and complete redistribution. The solution is performed using an iterative procedure, where monochromatic 
radiative transfer problems are successively solved. The present implementation is applicable for arbitrary model configurations 
with an optical depth up to 10 3 ~ 4 . Results of Lya line transfer calculations for a spherically symmetric model, a disk-like 
configuration, and a halo containing three source regions are discussed. We find the characteristic double-peaked Lya line 
profile for all models with an optical depth J; 1. In general, the blue peak of the profile is enhanced for models with infall 
motion and the red peak for models with outflow motion. Both velocity fields produce a triangular shape in the two-dimensional 
Lya spectra, whereas rotation creates a shear pattern. Frequency-resolved Lya images may help to find the number and position 
of multiple Lya sources located in a single halo. A qualitative comparison with observations of extended Lya halos associated 
with high redshift galaxies shows that even models with lower hydrogen column densities than required from profile fitting 
yield results which reproduce many features in the observed line profiles and two-dimensional spectra. 
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1. Introduction 

Hydrogen Lya as a prominent emission line of high redshift 
galaxies is important for the understanding of galaxy forma- 
tion and evolution in the early universe. Recent narrow-band 
imaging and spectroscopic surveys used the Lya line to iden- 
tify galaxies at very high redshift (e.g. Hu et al. 1998| Kudritzki 



ture of the objects and the fact that Lya is a resonance line 
require detailed radiative transfer modeling. 

The transfer of resonance line photons is profoundly de- 
termined by scattering in space and frequency. Analytical 



(Neufeld 1990 1 as well as early numerical methods (Adams 



1972; Hummer & Kunasz 1980) were restricted to one 



et al. p000| ; Rhoads et al. |2000| Fynbo et al. |2000[ [200 1| ). 

But beside from being a good redshift indicator, Lya emis- 
sion also bears information on the distribution and kinematics 
of the interstellar gas as well as the nature of the energy source. 
Lya observations of high redshift radio gal axies ( e.g. H ippelein 



dimensional, static media. Only recently, codes based on the 
Monte Carlo method were developed which are capable to in- 
vestigate t he more gene ral case of a multi-dimensional medium 
(Ahn etal.|200li^002). 



199rj , |!997t Villar- 



& Meisenhei mer 1993 ; van Ojik et al. 
Martin et al. 1999 ) reveal extended Lya halos with sizes > 
100 kpc which are aligned with the radio jet. The line profiles 
are often double-peaked and the two-dimensional spectra point 
to complex kinematics involving velocities > 1000 km s _1 . 

The interpretation of Lya observations is difficult, because 
high-redshift radio galaxies tend to be in the center of proto 
clusters, where the radio jet interacts with a clumpy environ- 
ment influenced by merging processes (e.g. Bicknell et al. 



In this paper, we introduce a finite element method 
for the solution of the resonance line transfer problem in 
three-dimensional, moving media and present the results of 
some simple, slightly optically thick model configurations. 
The basic, monochromatic code was originally developed by 
Kanschat (1996) and is described in Richling et al. (2001 



200C; Kurk et al. 2001 ). Actually, the three-dimensional struc- 
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hereafter Paper I). The three-dimensional method is particu- 
larly useful for scattering dominated problems in inhomoge- 
neous media. Steep gradients are resolved by means of an adap- 
tively refined grid. Here, we only specify the extension from 
the monochromatic to the polychromatic problem including the 
implementation of the Doppler-effect and complete redistribu- 
tion. 
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In Sect. 2, we review the equations for resonance line trans- 
fer in moving media. In Sect. 3, we describe some details re- 
garding frequency discretization and the form of the result- 
ing matrices for coherent scattering and complete redistribu- 
tion and give a short outline of the complete finite element al- 
gorithm. Then, the results of a spherically symmetric model 
(Sect. 4), a disk-like configuration (Sect. 5) and a model with 
three separate source regions (Sect. 6) are presented. A sum- 
mary is given in Sect. 7. 

2. Line Transfer in Moving Media 

We calculate the frequency-dependent radiation field in moving 
media by solving the non-relativistic radiative transfer equation 
in the comoving frame for a three-dimensional domain ft which 
in operator form simply reads 



(T + V + S + X (x, i/)) J(x, n, v) = /(x, v) 



(1) 



The transfer operator T, the "Doppler" operator T>, and the 
scattering operator S are defined by 



TT(x, n, v) 
2?Z(x, n, v) 

5Z"(x, n, v) 



n ■ W x X(k, n, v), 
d 

w(x, n) v—X(x, n, v) 
*(*) 

4.7T 



/ / i?(n, 0; n, v)X(x, n, v) doj dv. 
Jo Js 2 



Considering three dimensions in space, the relativistic invari- 
ant specific intensity X is six-dimensional and depends on the 
space variable x, the direction n (pointing to the solid angle duo 
of the unit sphere S 2 ), and the frequency v. 

The extinction coefficient x(x, v) = k(x, v) + u(x, v) is 
the sum of the absorption coefficient k(x, v) — n{x)(f)(y) and 
the scattering coefficient <r(x, v) = cr(x)0(i/). The frequency- 
dependence is given by a normalized profile function <f>. 
Usually, we adopt a Doppler profile 



i 



^/ttAv d 



exp 



Av D 



(2) 



where vq is the frequency of the line center. The Doppler width 
Avp and the Doppler velocity vd are determined by a thermal 
velocity fthcrm as well as a turbulent velocity wturb 



Av D 



-v D 



c 



therm 



^turb- 



(3) 



c is the speed of light. For situations, where wturb 3> Wtherm, the 
Doppler core is very broad and would dominate the Lorentzian 
wings of a Voigt profile. Then, the Doppler profile is a reason- 
able description of a line profile. Throughout this paper we use 
v D = l(T 3 c - 300 kins" 1 (see also Table [l]). 
For the source term 

/(x, v) - /fi(x, i/)B(T(x), v) + e(x, v), (4) 

we can consider thermal radiation and non-thermal radiation. 
In the case of thermal radiation, / is calculated from a temper- 
ature distribution T(x), where B(T, v) is the Planck function. 



The "Doppler" operator T> is responsible for the Doppler 
shift of the photons. A derivation of the operator for non- 
relativistic velocities (v/c < 0.1) can be found in Wehrse et 
al. ( 200C| ). In contrast to the full relativistic transfer equation 
(cf. Mihalas & Weibel-Mihalas 1984), we neglect any terms 
responsible for aberration or advection effects. The function 



-n ■ V T 



(5) 



is the gradient of the velocity field v(x) in direction n. Note 
that the sign of w may change depending on the complexity of 
the velocity field v. 

The scattering operator S depends on the general redistri- 
bution function R(n, v\ n, v) giving the probability that a pho- 
ton (n, v) is absorbed and a photon (n, v) is emitted. In the 
following, we assume isotropic scattering 



51= - 



o(x) 

47T 



J R(pi v) J X(x, n, 0) dujdv, (6) 



where R(v, v) is the angle-averaged redistribution function 

R(0, v) — .} .„ { { R(x,h 7 i>;n,v)dujduj. (7) 
(4tt)^ J S 2 J S 2 

The function defined by (M) is normalized such that 



v) dv dv = 1. 



(8) 



For the calculations in Sect. 3, we consider two limiting cases: 
strict coherence and complete redistribution in the comoving 
frame. In the former case, we have 



and in the latter 



R(P, v) = cj>(v)8{v-v) 



R(v, v) = <j>{v)<j>{v) 



(9) 



(10) 



Thus, for coherent isotropic scattering, the scattering term sim- 
plifies to 

5 coh X(x, n, v) = - a ^ ^ / X{^h 1 v)duj (11) 

47T J S 2 

In the case of complete redistribution, the photons are scattered 
isotropically in angle, but are randomly redistributed over the 
line profile. Then, the scattering term reads 

5 crd I(x, n,v)=- Y ^ / (f>(v) f Z(x, n, v) dcu dv. 
47r Jo Js 2 

(12) 
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3. The Finite Element Method 



and 



3.1. Boundary Conditions 

For the modeling of prominent resonance lines, in particular 
Lya, we restrict the frequency discretization to a finite interval 
A := [i/q, vn+i], where v$ and vn+i are located far out of the 
line center in the continuum. The function io(x, n) from Eq. 
(|5]) defines the Doppler shift of the spectrum towards smaller 
(w < 0) or larger (w > 0) frequencies in the comoving frame. 
Therefore, boundary conditions of the form 

l(x, n, v) = I con t(x, n, v) 

are necessary on the upper and lower frequency interval bound- 
ary of the whole computational domain 

£ = £1 x S 2 x {uq, v N+ i}. 

Furthermore, to be able to solve Eq. ([!]), boundary conditions 
of the form 

J(x,n,i/) =2: in (x, n, v) 

must be imposed on the "inflow boundary" 

FxA= {(x,n,z/) G r|n r ■ n < 0}, 

where nr is the unit vector perpendicular to the boundary sur- 
face r of the spatial domain f2. The sign of the product nr • n 
describes the "flow direction" of the photons across the bound- 
ary. If we neglect any continuum emission (X con t = 0) and 
assume that there are no light sources outside the modeled do- 
main as in the case of a non-illuminated atmosphere (T m = 0), 
the two boundary conditions for the solution of the transfer 
equation ([j]) in moving media are 



I(x, n, v) = on E, 
X(x, n, v) = on T~ x A. 

3.2. Discretization for Coherent Scattering 



(13) 
(14) 



In order to solve the radiative transfer equation ([j]), we first 
carry out the frequency discretization including coherent scat- 
tering. With the abbreviation 



A coh = J 

Eq. ([[]) can be written as 



^coh 



A coh l(x, n, v) 



w(x, n)v—l(x,n,v) 



(15) 



f(x,v). (16) 



In Paper I, we described a Galerkin discretization for the 
monochromatic radiative transfer equation. Considering mem- 
ory and CPU requirements this approach is by far too "expen- 
sive" for the frequency-dependent transfer problem. Instead, 
we use an implicit Euler method for N equidistantly distributed 
frequency points z/; 6 {vi, z^ 2 , vn} C A. Since the function 
w(x, n) may change its sign, this simple difference scheme for 
the Doppler term reads 



io(x, n)v 



81 
dv 



(wi > 0) 



(17) 



di 

w(x, n)v— 
ov 



1 Av 



{Wi < 0), 



(18) 



where Az/ is the constant frequency step size. All quantities 
referring to the discrete frequency point z/j are denoted by an 
index "i". Employing the Euler method, we get a semi-discrete 
representation of Eq. ( |l6| ) 



A 9 



\W\Vj 



^ = h + 



(w > 0) 
(w < 0) 



(19) 



The additional term on the left side of Eq. ( |19| ) can be inter- 
preted as an artificial opacity, which is advantageous for the 
solution of the resulting linear system of equations. The addi- 
tional term on the right side of Eq. ( |T^ ) is included as an arti- 
ficial source term. Equation ( |T9| ) can be written in a compact 
operator form 



A? h ^ = fi 



(20) 



Given a discretization with L degrees of freedom in il, M di- 
rections on the unit sphere S 2 and N frequency points, the 
overall discrete system has the matrix form 



A coh u = f , 



(21) 



with the vector u containing the discrete intensities and the 
vector f the values of the source term for all frequency points 
Vi. Both vectors are of length (L ■ M ■ N) and A coh is a 
(L ■ M ■ N) x (L ■ M ■ N) matrix. The fact that the func- 
tion w;(x, n) may change its sign, results in a block-tridiagonal 
structure of A coh and we get 



coll 



B 2 


V 



Ri 

A^ oh R 2 



\ 



'• Rjv-i 
B^v A™ h 



/ U! \ 
U 2 



\u N J 



/ f l\ 

f 2 



Vw 



According to the sign of w(x, n) the block matrices R,; and B,; 
hold entries of w(x, n)vi/Av causing a redshift and blueshift 
of the photons in the medium, respectively. Requiring a reason- 
able resolution, the resulting linear system of equations of the 
total system is too large to be solved directly. Hence, we are 
treating N "monochromatic" radiative transfer problems 



. coh. 



with a slightly modified right hand side 

f i = f i + RiUi + i + BiUi_i. 



(22) 



(23) 



Note that using simple velocity fields (e.g. infall or outflow) 
the sign of w(x, n) is fixed and the matrix A coh has a block- 
bidiagonal structure. This is favorable for our solution strategy, 
since we only need to solve Eq. (E2I) once for each frequency. 
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3.3. Discretization for Complete Redistribution 

Considering complete redistribution, we use an implicit Euler 
scheme to discretize the Doppler term as described above and 
a simple quadrature method for the frequency integral in the 
scattering operator <S crd in Eq. (|l2|). Starting from N equidis- 
tantly distributed frequency points Vi 6 \y\i •••> ^n} C A 
and N weights q-y, q^i ■■■> In, we define a quadrature method 



N 



for integrals f A £(v')dv'. In the case of complete redistribution 
the kernel is 



4-7T 



/ X(x, n, Vj)duj. (25) 
Js 2 



Separating the terms with the unknown intensities Xj from 
the known quantities Xj the discretized scattering integral ([IJ 
reads 



47T 



f cr- N f 

qi / lidui + } 4>jqj / I(x, n, Vj)du. (26) 
Js 2 47r ^ Js 2 



Using this quadrature scheme and employing the Euler method 
for the frequency derivatives we get a semi-discrete formula- 
tion of the transfer problem for each frequency point including 
complete redistribution 



A 



crd 



Au 



N 



— fi + — ^2 <f>jQj I > '- 1 

3& S 



T(x.,h,Vj)du), (27) 



where 



A? d = r + Xi + fas 



coh 



(28) 



The additional terms on the right hand side of Eq. ( 27 ) must 
be interpreted as artificial source terms. Eq. (27) can also be 
written in a compact operator form 



•A-i It ft- 



(29) 



The total discrete system has the matrix form (cf. Eq. (|21J)) 

f. (30) 



A crd u 



Unfortunately, the global frequency coupling via the scattering 
integral ( |l2| ) results in a full block matrix and we get 



. crd 



A crd 



R-i + Q2 Q.3 
B 2 + Qi Af d R 2 + Q 3 

Qi 



V Qi 



Qn \ 



A crd / 



B; and R s ; again contain the contribution of the Doppler factor 
w(x, r\)vi/&v, whereas Qj holds the terms from the quadra- 
ture scheme. As we already explained in the case of coherent 



scattering, we do not solve the total system, but N "monochro- 
matic" radiative transfer problems 

A crd _ i 



(31) 



with a modified right hand side 



fj = fj + RiU i+ i + BiUj-i + ^ Qj u j- ( 32 ) 



(24) 3.4. Full Solution Algorithm 



Eq. ( pop and (|29|) are of the same form as the monochromatic 
radiative transfer equation, AX — /, for which we proposed a 
solution method based on a finite element technique in Paper I. 
The finite element method employs unstructured grids which 
are adaptively refined by means of an a posteriori error estima- 
tion strategy. Now, we apply this method to Eq. (20) or Eq. ([29|). 
In brief, the full solution algorithm reads: 

1. Start with 1 = for all frequencies. 

2. Solve Eq. @ or Eq. @ for i = 1, .., AT. 

3. Repeat step 2 until convergence is reached. 

4. Refine grid and repeat step 2 and 3. 

We start with a relatively coarse grid, where only the most im- 
portant structures are pre-refined, and assure that the frequency 
interval \v\ , vn] is wide enough to cover the total line profile. 
Then, we solve Eq. ([20]) or Eq. ( p9| ) for each frequency sev- 
eral times depending on the choice of the redistribution func- 
tion and the velocity field. During this fix point iteration, we 
monitor the changes of the resulting line profile in a particular 
direction n out . Not until the line profile remains unchanged, 
we go over to step 4 and refine the spatial grid. Again, we ap- 
ply the fixed fraction grid refinement strategy: The cells are 
ordered according to the size of the local refinement indicator 
r\K — rnax{r\K(yi))\ Vi and a fixed portion of the cells with 
largest r\K is refined. r\K(yi) is an indicator for the error of the 
solution in cell K at frequency vi. 

In Paper I (Sect. 2.3), we introduced various possibilities to 
determine t]k- If we use the total flux leaving the computational 
domain in direction n out as the quantity whose error functional 
is used to calculate t\k (dual method), our convergence cri- 
terion described above is perfectly reasonable. We found that 
this criterion is also useful for the much simpler method (L 2 
method), where t\k is determined from the global error func- 
tional. In this case, the rate of convergence of the line profile 
in the directions n 7^ n out and for the line profile in direction 
n out are very similar. We are aware that this result may not be 
valid for more complex model configuration. 

4. Model Configurations 

4. 1. Halo and Sources 

We investigate spherically symmetric and elliptical distribu- 
tions of the extinction coefficient x(x) = x(x,y,z) of the form 

,.2\ 



Xo/(l 

X(x) = { X o/(l + ar 2 ) 

Xo /(l + ar2)/10 3 



for r < r c 
for r c < r < rh , (33) 
for r > rh 
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where r 2 = (x/a) 2 + (y/b) 2 + (z/c) 2 . Note, that the value of 
the extinction coefficient is constant in the center within core 
radius r c and outside halo radius r^. If not stated otherwise, xo 
is determined from the line center optical depth 

rn, 

xW^^o) n thick ax (34) 

between r c and ru along the most optically thick direction 
nthick- In total, the spatial distribution of x is determined by 
six parameters: the length of the half-axes a, b and c, the radii 
r c and 77,, the dimensionless parameter a, and the optical depth 
t. For r c , Th and a we use the values given in Table [j]. 

Since we are predominantly interested in the transfer of ra- 
diation in resonance lines like Lya, we assume cr(x) = x( x ) 
and «(x) = for all calculation presented here. This restricts 
us to the use of purely non-thermal source functions. In particu- 
lar, we consider one or several spatially confined source regions 
with radius r s each centered at a position x^: 

/(x,i/) = /JM !° r ! x - x <|^ . (35) 
J v ; \0 for |x - Xj| > r s 

The function 4>{v) is the Doppler profile defined in Eq. (|^). 

4.2. Velocity Fields 

In general, the finite element code is able to consider arbitrary 
velocity fields. For velocity fields defined on a discrete grid, for 
example for velocity fields resulting from hydrodynamical sim- 
ulations, the velocity gradient in direction n must be obtained 
numerically. Here, we use two simple velocity fields and cal- 
culate the function w analytically. 

The first velocity field describes a spherically symmetric 
inflow (vq < 0) or outflow (vq > 0) and is of the form 

v io = v (—) -, (36) 
V r / r 

where r = |x| and vq the scalar velocity at radius vq. The cor- 
responding w function is 



Table 1. Parameters used for all calculations. Distances are 



,(x,n)=, pyfi-(l + l)H 
\ r / \r r 



(37) 



For the second velocity field, we assume rotation around 
the z-axis 



v ro t = w 1 ^j^J R 1 | -x 



(38) 



where R 2 = x 2 + y 2 is the distance from the rotational axis 
and v the scalar velocity at distance R . If n = (n x ,n y , n z ), 
the w function reads 

w = V0 {-r) + W )■ 

(39) 

In the special case of a Keplerian velocity field (I = 0.5), we 
obtain 

w = v Q ^-R° 5 R~ 3 - 5 [xy(n 2 y - n\) + n x n y (x 2 - y 2 )} (40) 



(see also Baschek & Wehrse |1999[ ). Again, the values used for 
vq, ro and Rq are given in Table 111 



given in units of the computational cube (see Sect. 4.3 ) 





a r s v D 


v 


ro 


Ro 


1.0 0.2 


10 3 0.2 10~ 3 c 


-10" 3 c 


0.2 


1.0 



4.3. Normalization 

Since we use a normalized form of the equation of transfer, 
where the computational domain is the unit cube, the results 
are valid for model configurations with different length scales 
and for different resonance lines. When 2x max is the size of the 
unit cube, coordinates in physical units are simply obtained by 
the transformation x — > x • x max . 

A particular line and the corresponding frequency scale 
must be defined when transforming to the observers frame. 
From the solution X(x, n, v) in the comoving frame, we obtain 
the solution X(x, n, v) in the observers frame by performing 
the transformation 



-n 3 

2(x, n, v) = X(x, n, v) ( — 



where 



v 1 



v(x) 



(41) 



(42) 



The relation between optical depth r and column density 
of the scattering media depends on the central line absorption 
cross section x{ v o) an d therefore on the Doppler velocity vq. 
For the Lya line of neutral hydrogen the central line absorption 
cross section is 



XLya 



2.0 x 10" 



300 km s" 



cm 



(43) 



Hence, the column density of neutral hydrogen for a halo with 
given t is 



iV H 



0.5 x 10 14 t 



(sOOim^) Cm " (44) 



We would like to emphasize that this value is the column den- 
sity of the neutral halo alone. Considering the column density 
of the source region, the total column density is about twice 
as high and is in the range A^h ~ [10 13 , 10 16 ] cm -2 for the 
investigated configurations with r = [0.1, 100]. 

In comparison to column densities deduced from Voigt pro- 
file fitting procedures of Lya profiles of high redshift radio 
galaxies (van Ojik et al. 1997), our values are at least an or- 
der of magnitude too low. With the present implementation of 
our finite element code, we could calculate systems with col- 
umn densities up to Ah ~ 10 17 — 10 18 cm~ 2 , but only with a 
large amount of computational time. Still higher column densi- 
ties are beyond feasibility. 

5. Test calculations 

We start with a spherically symmetric model configuration 
(a = b = c = 1) and a single source region located at x = 0. 
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Fig. 1. Lya line profiles calculated with the FE code for a spherically symmetric model configuration: a) a static halo with 
coherent scattering, b) an infalling halo with coherent scattering, c) a static halo with complete redistribution, and d) an infalling 
halo with complete redistribution. For the static cases a) and c) the line styles refer to calculations with different optical depth r 
as indicated. The small window in a) enlarges the peak of the line. The crosses mark the results of the analytical solution. For the 
moving halos we show in b) the results for r = 1 (thin lines) and r = 10 (thick lines) and in d) only for r = 10. Here, the line 
styles refer to the exponent I used for the velocity fields. 



Fig. |l| shows the results of the finite element code (FE code) 
for different optical depths, velocity fields and redistribution 
functions. We used 41 frequencies equally spaced in the inter- 
val (y — vq)I A^d = [—4, 6] and 80 directions. Starting with a 
grid of 4 3 cells and a pre-refined source region, we needed 3-5 
spatial refinement steps. 

The simplest case is a static model with coherent isotropic 
scattering. Fig. |l]a displays the emergent line profiles for dif- 
ferent r. As expected, the Doppler profile is preserved and the 
flux F v is independent on r. The deviation of the numerical re- 
sults from the analytical solution indicated with crosses is very 
small. The line profiles for r = 0.1 and r = 1 are identi- 
cal even in the little window which shows the peak of the line 
in more detail. For r = 100, the total flux is still conserved 
better than 99%. This result demonstrates that the frequency- 
dependent version of our FE code works correctly. 

Next, we consider an infalling halo with coherent scattering 
and show the effects of frequency coupling due to the Doppler 
term. The emergent line profiles in Fig. [j]b are plotted for differ- 
ent exponents I of the velocity field Vi Q defined in Eq. (|36[). The 



line profiles are redshifted for r = 1 (thin lines). Most of the 
photons directly travel through the halo moving away from the 
observer. Since the Doppler term is proportional to the gradient 
of the velocity field, the redshift is larger for a greater exponent 
I. For t — 10 (thick lines), the line profiles are blueshifted. 
Before photons escape from the optically thick halo in front 
of the source, they are back-scattered and blueshifted in the 
approaching halo behind the source. The blueshift is less pro- 
nounced for the accelerated infall with I = 2, because the 
strong gradient of the velocity field leads to a slight redshift 
in the very inner parts of the halo. In this region, the total opti- 
cal depth is still small. Further out, where the total optical depth 
increases, the line profile becomes blueshifted. 

Complete redistribution is another method of frequency 
coupling which leads to a stronger coupling than the Doppler 
effect (see Sect. 3). The line profiles obtained for a static model 
with complete redistribution are displayed in Fig. []]; for differ- 
ent t. With increasing optical depth the photons more and more 
escape through the line wings. For r > 1, we get a double- 
peaked line profile with an absorption trough in the line center. 
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a) 



b) 
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Fig. 2. Results obtained for a static, disk-like model configuration: a) Emergent line profiles for different viewing angles and 
optical depths r = r(iithick)- b) Two-dimensional spectra for an edge-on view (90°) and a nearly face-on view (20°) and 
different optical depths. The position and width of the slit is indicated in Fig. |^a. The contours are given for 2.5%, 5%, 7.5%, 
10%, 20%, 90% of the maximum value. 



The greater r the larger the distance between the peaks and the 
depth of the absorption trough. Since our frequency resolution 
is too poor for the pointed wings, the flux conservation is only 
96% for r = 100. 

Fig. [j]d shows the results for an infalling halo with com- 
plete redistribution for r = 10 and different exponents I. For 
I = and I = 0.5 the infalling motion of the halo enhances the 
blue wing of the double-peaked line profile. Equally, an out- 
flowing halo would enhance the red peak. But for I = 2, the 
red peak is slightly enhanced for an infalling halo due to the 
strong velocity gradient, as explained above. This example af- 
fords an insight into the mechanisms of resonance line forma- 
tion and shows the necessity of a profound multi-dimensional 
treatment. 

6. Applications 

All calculations discussed in this section were performed with 
complete redistribution. We used 49 equidistant frequencies 
covering the interval (v — v a )//S.V£, = [—6,6], 80 directions 



and started from a grid with 4 3 cells and pre-refined source re- 
gions, which results in several 10 3 cells for the initial mesh. 
3-7 refinement steps were performed leading to approximately 
10 5 cells for the finest grid. 

6.1. Elliptical Halos 

As a first step towards a full three-dimensional problem without 
any symmetries, we investigated an axially symmetric, disk- 
like model configuration (a : b : c = 3 : 3 : 1) with a single 
source located at x = for several r = r(nthick) and for dif- 
ferent velocity fields. The most optical thick direction nthick is 
the direction within the equatorial plane of the disk. The direc- 
tion perpendicular to the equatorial plane is the z-axis which 
we also call rotational axis even for cases without rotation. 

By means of the static case we show which kind of infor- 
mation can be obtained from the results of our FE code. The 
emergent line profiles provide global information on the under- 
lying system. They are displayed in Fig. ||a for different optical 
depths and viewing angles. The viewing angle is defined as the 
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Fig. 3. Spatial intensity distribution for a static disk-like model configuration with T(n t hi c k) = 10 for different viewing angles: a) 
Frequency integrated intensity distribution, b) Intensity distribution for specific frequencies. The frequency is given in Doppler 
units relative to the line center (y — v$)/ A^d at the right hand side. 



angle between the rotational axis and the direction towards the 
observer. For instance, a viewing angle of 90° means, the disk 
is "observed" edge-on. Again, we obtain the characteristical 
double-peaked line profile for r <; 1. The optical thickness de- 
creases for smaller viewing angles. For a viewing angle of 0° 
the effective optical depth is only r(n t hick)/3. Consequently, 
the distance of the peaks in the line profile and the depth of the 
absorption feature is growing with increasing viewing angle. 



Note, that the total relative flux F(n = 0°)/F(n = 90°) es- 
caping along the z-axis is increasing with increasing r(n t hick)- 
Fig. ||b will be explained later. 

Far more information is contained in the spatial distribu- 
tion of the line intensity projected into the s-i-plane which is 
the plane perpendicular to the viewing direction. In Fig. ||a the 
frequency-integrated intensity distribution is shown for r = 10 
and different viewing angles. At 90° the Lya image shows an 
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Fig. 4. Results obtained for a disk-like model configuration with global infall motion (7 = 0.5): a) Emergent line profiles for 
different viewing angles and optical depths r = r(n t hick)- b) Two-dimensional spectra for an edge-on view (90°) and a nearly 
face-on view (20°) and different optical depths. The position and width of the slit is indicated in Fig. [|a. The contours are given 
for 2.5%, 5%, 7.5%, 10%, 20%, 90% of the maximum value. 



elliptical emission line region. The diffuse halo of scattered ra- 
diation is most clearly visible. With decreasing viewing angle 
the diffuse halo becomes less pronounced. For a nearly face- 
on view (20°), the halo is optically thinner and the intensity 
distribution appears spherically symmetric. 

Fig. ||b displays the spatial distribution of the Lya inten- 
sity for different viewing angles (columns) at selected frequen- 
cies (rows). For the edge-on view, the intensity distribution 
strongly depends on the frequency. In the outer line wings at 
(y — vq)//S.ve, = ±2 the contours are elliptical and repro- 
duce the disk-like form of the source region. The peaks of the 
line profile are at (y — vq)/A.vd ~ ±1. Here, the major axis 
of the elliptical intensity distribution in the inner parts of the 
model is parallel to the rotational axis, whereas the major axis 
of the elliptical contour lines in the outer parts remains paral- 
lel to the equatorial plane of the disk. In the very center of the 
line ({v — Uq)/ A^d = 0), the direct view to the source region 
is blocked by optically thick material. There appear two knots 
of Lya emission directly above and below the disk, which are 
surrounded by an extended elliptical halo. The Lya knot below 



the disk vanishes with decreasing viewing angle. At 20° the 
intensity distribution is spherically symmetric for all frequen- 
cies. The observable emission region is most extended at the 
frequency of the line center. 

Two-dimensional spectra from high-resolution spec- 
troscopy provide frequency-dependent data for only one spatial 
direction. To be able to compare our results with these obser- 
vations we calculated two-dimensional spectra using the data 
within the slits indicated in Fig. |]a for the edge-on and nearly 
face-on view. The results are shown in Fig. ^b for different op- 
tical depths. We find that the form of the two-dimensional spec- 
tra is relatively insensible to the width of the slit. For r = 1 a 
single emission region is visible. But already at 90°, the highest 
contour lines reproduce the faint absorption trough of the cor- 
responding line profile (Fig. |^). The higher the optical depth 
the wider the gap and the spatial extent of the two peaks. Note 
that the spatial extent of the outer contour line only depends on 
T(rithick) and not on the viewing direction. 

What changes when imposing a macroscopic velocity field 
is shown in the following two figures. First, we consider an in- 
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Fig. 5. Results obtained for a disk-like model configuration with Keplerian rotation (/ = 0.5) around the z-axis: a) Emergent line 
profiles for different viewing angles and optical depths t = r(rithick). b) Two-dimensional spectra for an edge-on view (90°) and 
a nearly face-on view (20°) and different optical depths. The position and width of the slit is indicated in Fig. ||a. The contours 
are given for 2.5%, 5%, 7.5%, 10%, 20%, 90% of the maximum value. 



falling velocity field with / = 0.5 suitable for a gravitational 
collapse. Fig. ||a displays the calculated line profiles for dif- 
ferent r and viewing angles. As expected, the blue peak of 
the line is enhanced. The higher r the stronger the blue peak. 
Fig. [|b shows the corresponding two-dimensional spectra ob- 
tained with the same slit width and position as in the static case. 
For low optical depth, the shape of the contour lines is a tri- 
angle. Photons changing frequency in order to escape via the 
blue wing are also scattered in space. The consequence is the 
greater spatial extent of the blue wing. Apart from a growing 
gap between the two peaks with higher optical depth, the gen- 
eral triangular shape in conserved. 

Next, we investigated the elliptical model configuration 
with Keplerian rotation (/ = 0.5), where the z-axis is the axis 
of rotation. The results are plotted in Fig. || for r = 1 and 
r = 10. The emergent line profiles are symmetric with respect 
to the line center and show the same behavior with growing op- 
tical depth as in the static case. However, the extension of the 
line wings towards higher and lower frequencies is strongly in- 
creasing with increasing viewing angle because of the growing 
effect of the velocity field. Rotation is clearly visible in the two- 
dimensional spectra (Fig. ^p). The shear of the contour lines 
is the characteristical pattern indicating rotational motion. For 
an edge-on view at r = 10, Keplerian rotation produces two 
banana-shaped emission regions. 

In spite of the relatively low optical depth of the simple 
model configurations, our results reflect the form of line pro- 
files and the patterns in two-dimensional spectra of many high 
redshift galaxies. For example, the two-dimensional spectra 



of the Lya blobs discovered by Steidel et al. (200C, Fig. 8) 
are comparable to the results obtained for infalling (Fig. [|b) 
and rotating (Fig. ^p) halos. In the sample of van Ojik et 
al. (1997) are many high redshift radio galaxies with single- 



peaked and double-peaked Lya profiles. The corresponding 
two-dimensional spectra show asymmetrical emission regions 
whi ch ar e more or less extended in space. De Breuck et 
al. ( |2000| ) find in their statistical study of emission lines from 
high redshift radio galaxies that the triangular shape of the Lya 
emission is a characteristical pattern in the two-dimensional 
spectra of high redshift radio galaxies. Since the emission of the 
blue peak of the line profile is predominately less pronounced, 
most of the associated halos should be in the state of expansion. 



6.2. Multiple Sources 

High redshift radio galaxies are found in the center of proto 
clusters. In such an environment, it could be possible that the 
Lya emission of several galaxies is scattered in a common 
halo. To investigate this scenario, we started with a spheri- 
cally symmetric distribution for the extinction coefficient and 
with three source regions located at x\ = [0.5, 0.25, 0], X2 = 
[-0.5, 0.25, 0] and x 3 = [0, -0.25, 0] forming a triangle in the 
x-y-plane. In the following figures, the results for r = 10 are 
presented. 

We begin with the static model. Figure ^a shows Lya im- 
ages for four selected viewing directions. The specified angle 
is the angle between the viewing direction and the x-y-plane. 
Note that the orientation of the source positions within the 
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Fig. 6. Spatial intensity distribution for a static spherically symmetric model configuration with t — 10 and three source regions 
for different viewing angles: a) Frequency integrated intensity distribution, b) Intensity distribution for specific frequencies. The 
frequency is given in Doppler units relative to the line center (y — isq)/ A^d at the right hand side. 



plane is different for each image. Viewing the configuration 
almost perpendicular to the x-y-plane (70°), renders all three 
source regions visible, because the source regions are situated 
in the more optically thin, outer parts of the halo. Remember 
that most of the scattering matter is in and around the center 
of the system. For other angles, some of the source regions 
are located behind the center and therefore less visible on the 
images. Figure ||b shows images at different frequencies. In 



the line wings at [y — i/q)/Ai/£> = ±2, the number and po- 
sition of the source regions can be determined for all view- 
ing angles. However, at frequencies around the line center at 
(y — vq)I Ai>d = ±1 and 0, the number of visible sources 
depends on the viewing angle. Some of the images show only 
one, other two or three sources. 

The corresponding line profiles and two-dimensional spec- 
tra are displayed in Fig. [7] for the static case as well as for a 
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Fig. 7. Results obtained for a spherically symmetric model configuration with r = 10 and three source regions for the static case 
and two different velocity fields (inflow and rotation): a) Emergent line profiles for different viewing angles, b) Two-dimensional 
spectra for a nearly face-on view (70°) and a nearly edge-on view (10°). The position and width of the slit is indicated in Fig. |6^. 
The contours are given for 2.5%, 5%, 7.5%, 10%, 20%, 90% of the maximum value. 



halo with global inflow and Keplerian rotation. Width and po- 
sition of the slits are depicted in Fig. ||a. We get double-peaked 
line profiles for almost all cases, except for the rotating halo, 
where the line profiles are very broad for viewing angles lower 
than 70°. Additional features, dips or shoulders, are visible in 
the red or blue wing. They arise because the three sources have 
significantly different velocities components with respect to the 
observer. 



This example demonstrates the complexity of three dimen- 
sional problems. In a clumpy medium, the determination of 
the number and position of Lya sources would be difficult by 
means of frequency integrated images alone. Two-dimensional 
spectra may help, but could prove to be too complicated. More 
promising are images obtained from different parts of the line 
profile or information from other em ission lines of OIII or Ha 
in a manner proposed by Kurk et al. (2001). 



The slit for a viewing angle of 70° contains two sources. 
They show up as four emission regions in the two-dimensional 
spectra (Fig. ^|b). The pattern is very symmetric, even for the 
moving halos. For a viewing angle of 10°, the slit covers all 
source regions. Nevertheless, the two-dimensional spectra are 
dominated by two pairs of emission regions resulting from the 
sources located closer to the observer. The third source region 
only shows up as a faint emission in the blue part in the case of 
global inflow. In the case of rotation, emission regions from a 
third source are present. But the overall pattern is very irregular 
and prevents a clear identification of the emission regions. 



7. Summary 

We presented a finite element code for solving the resonance 
line transfer problem in moving media. Non-relativistic ve- 
locity fields and complete redistribution are considered. The 
code is applicable to any three-dimensional model configura- 
tion with optical depths up to 10 3-4 . 

We showed applications to the hydrogen Lya line of 
slightly optically thick model configurations (r < 10 2 ) and 
discussed the resulting line profiles, Lya images and two- 
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dimensional spectra. The systematic approach from very sim- 
ple to more complex models gave the following results: 

- An optical depth of r <; 1 leads to the characteristic dou- 
ble peaked line profile with a central absorption trough as 
expected from analytical studies (e.g. Neufeld 199C). This 
form of the profile is the result of scattering in space and 
frequency. Photons escape via the line wings where the op- 
tical depth is much lower. 

- Global velocity fields destroy the symmetry of the line pro- 
file. Generally, the blue peak of the profile is enhanced 
for models with infall motion and the red peak for models 
with outflow motion. But there are certain velocity fields 
(e.g. with steep gradients) and spatial distributions of the 
extinction coefficient, where the formation of a prominent 
peak is suppressed. 

- Double-peaked line-profiles show up as two emission re- 
gions in the two-dimensional spectra. Global infall or out- 
flow leads to an overall triangular shape of the emission. 
Rotation produces a shear pattern resulting in banana- 
shaped emission regions for optical depths <; 10. 

- For non-symmetrical model configurations, the optical 
depth varies with the line of sight. Thus, the total flux, the 
depth of the absorption trough and the pattern in the two- 
dimensional spectra strongly depend on the viewing direc- 
tion. 
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The applications demonstrate the capacity of the finite el- 
ement code and show that the three-dimensional structure and 
the kinematics of the model configurations are very important. 
Thus, beside exploring higher optical depths up to the limits 
of our code, we will consider clumpy density distributions as 
well as dust absorption and try to model the Lya emission of 
individual high redshift galaxies. In addition, we intend to im- 
plement a second order method for the frequency discretiza- 
tion. Furthermore, we plan to overcome the difficulties in solv- 
ing the line transfer problem for configurations with r > 10 4 . 
Therefore, we will extend our method and use the diffusion ap- 
proximation in the most optically thick regions. In the other re- 
gions, the full frequency-dependent line transfer problem must 
be solved. 
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